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Abstract 


The  theoretical  basis  for  the  Landau-Teller  equation 
commonly  used  to  model  rotational  nonequilbrium  is  re¬ 
viewed.  Several  assumptions  underlying  this  model  are  in¬ 
dicated  to  be  unrealistic  for  rarefied  hypersonic  flow.  A  new 
rotational  nonequilibrium  model  based  on  recent  measure¬ 
ments  up  to  2000 K  of  quantum  state  transition  rates  is  pre¬ 
sented.  The  new  model  is  applied  to  the  continuum  study  of 
normal  shock  wave  structure  using  both  the  Burnett  equa¬ 
tions  and  a  simplified  nonlinear  constitutive  relation.  Com¬ 
parisons  are  made  between  shock  wave  temperature  profiles 
generated  from  the  new  model  and  those  generated  using 
the  Landau-Teller  model.  Comparisons  of  shock  recipro¬ 
cal  thicknesses  between  experimental  data  for  nitrogen  and 
continuum  solutions  using  both  rotational  models  are  made. 
The  new  rotational  model  agrees  well  with  experiment  up  to 
Mach  6,  and  under  predicts  shock  thickness  at  higher  Mach 
numbers.  The  Landau-Teller  model  agrees  well  with  exper¬ 
imental  shock  thickness  up  to  Mach  3,  and  over  predicts  the 
thickness  at  higher  Mach  numbers.  A  modification  to  the 
rotational  collision  number  in  the  Landau-Teller  model  is 
found  to  give  results  which  agree  with  experimental  shock 
thicknesses  at  all  Mach  numbers  up  to  1 1. 

Nomenclature 


Ai,  Bit  C,  =curve  fit  constants 
AOTV  =  Aero-assisted  Orbital  Transfer 
Vehicle 

(V  =  rotational  specific  heat 

Cyf  =  translational  specific  heat 
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=  one-dimensional  substantial  derivative 
=  rotational  energy  per  unit  mass 
=  effective  rotational  energy  of  translation 
-  reference  rotational  energy 
=  translational  energy  per  unit  mass 
=  pCvtTt  +  Er  +  1/2  puZ  , 
total  fluid  energy  per  unit  volume 
=  per,  rotational  energy  per  unit  volume 
=  an  energy  difference  per  unit  mass 
=  energy  of  ith  quantum  state 
=  non-dimensional  energy  of  ith  state 
=  non-dimensional  function 
=  flux  vector 

=  non-dimensional  function 
=  non-dimensional  function 
=  maximum  allowable  quantum  number 
=  Boltzmann’s  constant 
=  transition  rate  from  state  « to  / 

=  temporal  index 

=  an  exponent  for  the  new  rotational  model 
=  National  Aerospace  Plane 
=  number  density  of  itk  quantum  state 
=  pRTt,  thermodynamic  pressure 
=  one  dimensional  monatomic  heat  flux  vector 
=  addition  to  heat  flux  for  diatomic  gas 
=  gas  constant 
=  shock  reciprocal  thickness 
=  normalized  temperature 
=  rotational  temperature 
=  viscosity  reference  temperature 
=  translational  temperature 
=  reference  temperature  used  for  Kl=>) 

=  reference  temperature  used  for  Zr 
=  temperature  ratio 
=  function  of  T 
=  fluid  velocity 
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conservation  vector 
source  term  vector 
rotational  energy  source  term 
rotational  collision  number  for 
Landau-Teller  model 
infinite  temperature  collision  number 
rotational  collision  number  for  new 
rotational  model 
parameters  for  Kl=>) 

Shock  density  temperature  separation 
/itv,  rotational  thermal  conductivity 
ypcvt.  translational  thermal  conductivity 
viscosity 
viscosity  at  Tref 
viscosity-temperature  exponent 
one  dimensional  stress  tensor 
density 

normalized  density 
characteristic  rotational  temperature 
ir/r/4  p,  mean  collision  time 

Introduction 

Simulating  the  flow  about  planned  transatmospheric  ve¬ 
hicles,  such  as  AOTV  and  NASP,  and  about  hypersonic 
missies,  requires  the  inclusion  of  rotational  thermal  non¬ 
equilibrium  in  rarefied  flow  codes.  This  is  due  to  the  shock 
wave  thickness  becoming  a  sizable  portion  of  the  flowfield 
at  sufficiently  high  altitudes.  Chapman,  et  al  have  shown 
that  the  Navier-Stokes  equations  require  the  inclusion  of  ro¬ 
tational  thermal  nonequilibrium  as  well  as  modified  consti¬ 
tutive  relations  to  achieve  experimentally  measured  shock 
thicknesses  in  diatomic  gases  such  as  Ni .  Furthermore, 
they  show  that  the  inclusion  of  rotational  thermal  nonequi¬ 
librium  is  equally  as  important  as  the  modifing  the  constitu¬ 
tive  relations;  whereas,  Fiscko  and  Chapman2  have  shown 
that  modifying  the  constitutive  relations  alone  is  sufficient 
to  obtain  good  agreement  with  measured  shock  thicknesses 
in  monatomic  gases. 

Computational  models  for  internal  energy  relaxation  are 
usually  based  on  the  the  well-known  Landau-Teller  equa¬ 
tion.  The  model  has  been  widely  applied  to  both  vibrational 
and  rotational  nonequilibrium.  The  Landau-Teller  formula¬ 
tion  in  terms  of  rotational  relaxation  is 

deT  _e?(Tt)-er 
dX  '  ZRrc  ■ 

In  continuum  hypersonic  CFD  this  equation  is  generally 
used  directly  as  the  computational  algorithm.  In  particu¬ 
late  CFD  (e.g.,  the  Direct  Simulation  Monte  Carlo  (DSMC) 


technique),  it  is  used  indirectly  by  handling  the  transfer  of 
internal  and  translational  energy  during  a  molecular  colli¬ 
sion  with  the  technique  outlined  by  Larsen  and  Borgnakke.3 
This  treatment  of  internal  energy  relaxation  is  consuicted 
to  give  results  which  conform  closely  to  the  Landau-Teller 
model.  Unfortunately,  as  discussed  below,  the  Landau- 
Teller  equation  is  derived  under  relatively  restrictive  as¬ 
sumptions  which  are  not  fully  realistic  for  hypersonic  flow 
conditions.  This  is  the  case  for  both  vibrational  and  rota¬ 
tional  energy  relaxation. 

The  principal  shortcomings  of  the  Landau-Teller  model 
when  applied  to  vibrational  relaxation  are  embodied  in  the 
underlying  assumptions  (1)  that  a  vibrating  diatom  can  be 
modeled  as  a  harmonic  oscillator  hence  allowing  quantum 
jumps  of  only  one  state  in  a  collision,  and  (2)  that  disso¬ 
ciation  does  not  affect  vibrational  relaxation.  These  as¬ 
sumptions  are  overly  restrictive  when  considering  tne  high 
temperature  conditions  in  hypersonic  flows.  In  attempts 
to  overcome  these  shortcomings,  models  have  been  devel¬ 
oped  which  account  for  more  “diffusion-like”  (i.e.,  allow¬ 
ing  for  quantum  jumps  greater  than  one)  relaxation  (e.g.. 
Keck4 )  and  for  removal  of  highly  excited  vibrational  states 
through  dissociation  of  the  diatomic  molecules  (e.g..  Mar- 
rone  and  Treanor5 ).  For  vibrational  relaxation,  therefore, 
departures  from  Landau-Teller  relaxation  have  long  been 
recognized,  and  alternative  models  have  been  used  in  some 
recent  hypersonic  flow  computations  (e  g.,  Candler  and 
MacCormack6). 

However,  for  the  case  of  rotational  thermal  nonequi¬ 
librium,  we  are  not  aware  of  any  computational  algorithm 
being  used  in  continuum  CFD  other  than  Landau-Teller. 
Wang  Chang  and  Uhlenbeck7  have  shown  that  the  Landau- 
Teller  equation  can  be  derived  for  an  internal  energy  mode 
out  of  equilibrium  with  translational  energy  if  the  following 
conditions  apply: 

(1)  Internal  ;  nd  translational  temperatures  are  nearly  equal. 

(2)  In  „'na)  r  translational  energy  exchange  is  very  slow. 

(3)  Inter  tergy  is  in  a  Boltzmann  distribution. 

For  rotational  energy,  assumption  (1)  is  unrealistic  for  rar¬ 
efied  hypersonic  flow,  and  assumption  (2)  is  unrealistic  for 
diatomic  gases  such  as  N2,  O2,  and  NO,  for  which  ro¬ 
tational/translational  energy  exchange  is  not  exceedingly 
slow.  Thus,  while  assumptions  underlying  the  derivation 
of  the  Landau-Teller  model  for  rotational  nonequilibrium 
are  of  questionable  accuracy  for  rarefied  hypersonic  flow 
in  diatomic  molecules,  a  validation  of  the  Landau-Teller 
model  for  such  conditions  has  not  been  presented,  nor  has 
a  more  realistic  alternative  been  advanced  for  use  in  hyper¬ 
sonic  CFD. 
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It  should  be  noted  that  existing  determinations  of  ZR,  the 
number  of  collisions  for  rotational  relaxation,  generally  as¬ 
sume  the  Landau-Teller  equation  to  apply.  Such  is  the  case 
for  the  theoretical  determination  of  ZR  by  Parker*  (Eq.  (2)). 

7W 

z  =  _ ±K _  (2) 

1  +  (T’*/rt)1/2  +  (f  +  *)  (r*/rt) 


Development  of  New  Rotational 
Relaxation  Model 

The  method  used  to  develop  the  new  algorithm  involves 
solving  the  master  equation  for  stationary  adiabatic  rota¬ 
tional  relaxation.  The  master  equation  (Eq.  (3))  is  a  set 
of  rate  equations  which  describe  the  evolution  of  rotational 
quantum  state  populations. 


Lordi  and  Mates9  also  make  the  same  assumption  in  their 
computation  of  ZR.  Table  1  summarizes  the  parameters 
supplied  for  Eq.  2  by  Parker  and  the  values  required  to 
fit  the  Parker  model  to  the  Lordi  and  Mates  computations. 
In  the  case  of  experimental  determinations,  some  fluid  dy¬ 
namic  measurement  such  as  sound  absorption/dispersion, 
shock- wave  profiles,  or  free-jet  expansion  profiles  are  used 
to  deduce  ZR  from  the  Landau-Teller  equation.  Further¬ 
more,  most  of  the  experimental  data  are  limited  to  relatively 
low  temperatures  compared  to  those  encountered  in  hyper¬ 
sonic  flight.  Fig.  1  summarizes  the  analytic,  computational, 
and  experimental  values  for  ZR.  The  data  are  widely  scat¬ 
tered  at  room  tempertature  and  do  not  extend  beyond  about 
1200  Kelvin. 

In  view  of  the  questionable  theoretical  basis  of  the 
Landau-Teller  model  and  the  relative  uncertainty  of  ZR  at 
even  low  temeratures,  it  appears  that  the  appropriate  macro¬ 
scopic  model  for  treating  rotational  thermal  nonequilibrium 
in  high  temperature  rarefied  flows  is  uncertain.  There¬ 
fore,  this  paper  proposes  a  new  approach  for  the  compu¬ 
tation  of  rotational  energy  relaxation.  It  is  founded  on 
some  relatively  recent  measurements  by  the  inverse  Ra¬ 
man  spectroscopy  method  of  rotational  quantum  state  tran¬ 
sition  probabilities  for  nitrogen.  Objectives  of  the  present 
research  are: 

(1)  To  develop  from  the  spectroscopic  data  a  new  macro¬ 
scopic  model  for  rotational  relaxation. 

(2)  To  employ  the  new  model  and  the  Landau-Teller  model 
in  a  continuum  simulation  of  shock  wave  structure  in 
nitrogen,  a  flow  with  a  large  amount  of  rotational  non- 
equilibrium.  Two  different  constitutive  relations  will 
be  used:  (1)  the  Burnett  equations  shown  by  Fiscko 
and  Chapman10  to  yield  accurate  shock  profiles  for 
monatomic  gases,  and  (2)  a  simplified  nonlinear  con¬ 
stitutive  relation  based  on  earlier  work  of  Chapman,  et 
al.,1 

(3)  To  evaluate  the  accuracy  of  the  both  rotational  models  by 
comparing  computed  shock  reciprocal  thicknesses  with 
experimentally  measured  reciprocal  thicknesses  for  ni¬ 
trogen. 


-  Yi  Ki+jNi  ;  *  =  1 , ....  /max  ( 3) 

Eq.  (4)  represents  an  empirically  fit  relation  presented  by 
Rahn,  ef  a/.,111213  for  the  upward  ( i  <  j)  transition  rates. 


Ki^j  =  apF(Tt) 


( 1+jBj/gN 

V  1  +  fi  ) 


where: 


(4a) 


1-5 Ej  _  1.5 0ri(»+  1) 
kTt  ~  Tt 
A£,;  =  Ei  -  Ej 


F(Tt )  = 


(  I-exp(-m)  \ 
-  exp  {-mTt/T0) ) 


(46) 


The  downward  transition  rates  ( »  >  ;)  are  obtained  from 
the  upward  rates  by  microscopic  reversibility  (Eq.  (5)). 


2;  +  1 
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(5) 


Table  2  gives  the  values  suggested  by  Rahn,  et  a/.,”1213 
for  use  with  Nj  in  Eq.  (4).  Solution  of  the  master  equa¬ 
tion  yields  the  population  densities  as  a  function  of  time, 
and  by  summing  over  quantum  states,  as  shown  in  Eq.  (6), 
rotational  energy  as  a  function  of  time  for  the  adiabatic  re¬ 
laxation  process  is  achieved. 


er  =  £:»(»+  l)kerMi 

i-0 


Since  the  process  is  adiabatic,  total  energy  is  constant,  and 
translational  energy  is  easily  obtained  by  subtracting  rota-  1  — 
tional  energy  from  the  prescribed  total  energy.  Tempera¬ 
tures  are  defined  in  terms  of  energies  as  shown  by  Eq.  (7). - 
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Fig.  2  compares  temperature-time  relaxation  from  master 
equation  and  solution  from  the  Landau-Teller  equation  us¬ 
ing  Parker’s  model  for  Zr.  Here  the  parameters  in  Parker’s 
model  are  adjusted  to  fit  the  Lordi  and  Mates  computations 
(see  Table  l). 

It  is  from  these  adiabatic  relaxation  simulations  that  a 
new  single  rate  equation  model  for  rotational  relaxation 
can  be  inferred.  A  new  macroscopic  model  for  rotational 
thermal  nonequilibrium  is  proposed  having  the  following 
form: 


deT  =  iser  e'r  /|Aer|\" 
&  Zr  rc  I A  eT\  \  e'r  ) 

where: 


(8a) 


2 

e r  =  j(et+  er) 

A er  =  e?(Tt)  -  er 


(8b) 


Note  that  the  reference  energy,  e£,  represents  the  value  of 
rotational  energy  that  would  be  present  if  the  total  thermal 
energy  present  were  equipartitioned.  This  choice  makes 
the  reference  energy  constant  for  the  adiabatic  relaxation 
process,  but  it  will  not  be  constant  in  an  actual  applica¬ 
tion  such  as  the  shock  wave  problem.  The  form  of  Eq.  (8) 
is  such  that  it  is  applicable  to  both  situations  when  TV  is 
less  than  TV  and  when  TV  is  greater  than  TV.  This  corre¬ 
sponds  to  both  compressing  and  expanding  flows.  It  is  as¬ 
sumed  a  priori  that  Zr  and  n  will  not  be  constant  but  in 
general  functions  of  both  the  Tt  and  TV.  The  conventional 
Landau-Teller  method  is  recovered  if  one  chooses  n  =  1 
and^H(Ti,rr)  =  Zr(  TV),  a  function  of  translational  tem¬ 
perature  only. 

The  parameters  Zr  and  n  were  determined  as  follows. 
As  mentioned  above,  the  master  equation  solution  yields 
the  time  rate  of  change  of  rotational  energy.  Since  the  solu¬ 
tion  procedure  used  was  a  numerical  technique  for  solving 
stiff  systems  of  ordinary  differential  equations,  one  obtains 
the  energies  at  a  discrete  number  of  points  in  time  rather 
than  continuously  over  time.  The  temporal  derivative  of  the 
rotational  energy  was  then  obtained  by  summing  the  mas¬ 
ter  equation  (Eq.  (3))  over  quantum  space  at  each  discrete 
point,  l,  in  time. 


In  order  to  obtain  the  exponent  in  Eq.  (8a)  values  of  the 
energy  derivative  at  both  l  and  l  —  1  are  used  to  eliminate 
Zr  by  subtracting  the  common  logarithms  of  Eq.  (8a)  eval¬ 
uated  at  the  above  mentioned  points.  The  resulting  equa¬ 
tion  is  then  solved  for  n  as  shown  by  Eq.  (10)  for  the  case 
Ae,  >  0. 
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This  assumes  that  both  nand  Zr  are  only  weak  functions  of 
Tt  and  TV.  This  assumption  will  be  justified  by  the  results. 
Now  that  n  is  known,  Zr  is  obtained  by  simple  substitution 
into  Eq  (8a). 

Zr  and  n  were  calculated  from  numerous  solutions  of 
the  master  equation  where  two  parameters,  total  energy  and 
the  ratio  TV /7V ,  were  varied.  Initial  rotational  quantum  state 
distributions  were  assumed  to  be  Boltzmann  at  the  temper¬ 
ature  7V.  The  values  of  ZR  and  n  were  calculated  only 
at  points  near  the  initial  condition  since  it  was  found  that 
as  the  rotational  energy  departed  from  Boltzmann  distri¬ 
butions  the  values  of  Zr  and  n  were  no  longer  properties 
of  macroscopic  quantities,  but  began  to  depend  on  the  dis¬ 
tribution  itself.  This  assumption  unfortunately  limits  the 
new  model;  however,  as  mentioned  above,  the  derivation 
of  Wang  Chang  and  Uhlenbeck  also  makes  this  assump¬ 
tion.  Furthermore,  ultrasonic  determinations  of  Zr  and  the 
spectroscopic  determinations  of  K^}  are  made  using  data 
that  is  taken  when  a  Boltzmann  distribution  over  rotational 
states  exists. 

In  order  to  apply  the  above  result*;  to  hypersonic  CFD, 
it  is  desirable  to  obtain  suitable  curve  fits  for  n  and  Zr. 
The  equations  below  present  curve  fits  used  in  the  shock 
simulations  presented  in  this  paper. 
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G(T)  =Ax  +  A2(\T\A'  -  1) 
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+  Ttft  j  +  1 

(116) 


Table  3  summarizes  the  numerical  values  Ai,Bi,  and  Ci 
required  to  fit  Eq.  (11)  to  n  and  Zr.  These  curve  fits  were 
arrived  at  by  simple  trial  and  error  and  are  probably  more 
complicated  than  required,  they  however  do  accurately  fit 
the  calculated  values  for  n  and  Zr  for  temperatures  below 
5000  Kelvin  and  for  a  wide  range  of  Tr/Tt,  including  both 
Tr  <  T,  and  Tr  >  Tt. 

Figs.  3,  4,  5,  and  6  show  both  calculated  and  curve  fit 
values  of  n  and  Zr  as  function  of  Tt  for  selected  values  of 
Tr/Tt.  The  exponent  n  departs  from  Landau-Teller  value 
(n=l)  for  high  temperatures  and  a  high  degree  of  thermal 
nonequilibrium  (i.e.,  \Tt  -  TT\/Tt  »  I).  The  exponent 
n  also  behaves  in  a  fundamentally  different  manner  for 
temperature  ratios  corresponding  to  expanding  flows  com¬ 
pared  with  temperature  ratios  corresponding  to  compress¬ 
ing  flows.  The  exponent  seems  to  asymptote  to  constant 
value  (with  respect  to  Tt)  of  1.128  as  the  temperature  ra¬ 
tio  Tr/Tt  approaches  the  limiting  values  of  one  from  both 
compressing  and  expanding  values.  The  fact  that  this  value 
is  only  1 3%  different  than  the  Landau-Teller  value  suggests 
that  perhaps  the  Landau-Teller  equation  may  be  sufficient  as 
a  rough  model  for  rotational  nonequilibrium  even  though 
the  theoretical  basis  of  the  model  is  questionable.  Also 
note  that  n  and  Zr  are  not  strongly  dependent  on  Tt  and 
Tr  justifying  the  assumption  necessary  to  calculate  n  and 
ZR  above. 

Constitutive  Relations  for  Use 
in  Shock  Wave  Structure 

It  has  long  been  known  that  the  Navier-Stokes  equations 
are  inaccurate  for  flows  that  contain  strong  departures  from 
translational  equilibrium.  Such  flows  are  characterized  by 
large  changes  in  flow  quantities  over  the  space  of  a  few 
mean  free  paths.  Normal  shock  structure,  where  flow  quan¬ 
tities  such  as  translational  temperature  can  jump  by  as  much 
as  two  orders  of  magnitude  in  the  space  of  ten  or  so  mean 
free  paths,  is  one  such  flow  where  the  Navier-Stokes  equa¬ 


tions  yield  erroneous  results. 

In  the  case  of  shock  wave  structure,  the  shortcomings 
of  the  Navier-Stokes  equations  are  the  constitutive  rela¬ 
tions,  not  the  continuum  description  of  the  fluid.  This 
has  been  shown  by  Fiscko  and  Chapman10  for  monatomic 
gases  by  utilizing  the  Burnett  equations  to  solve  for  shock 
wave  structure  and  comparing  the  results  with  shock  struc¬ 
tures  from  a  non-continuum  DSMC  calculation.  Fiscko 
and  Chapman  found  that  the  Burnett  equations  yield  sig¬ 
nificantly  improved  continuum  solutions  of  shock  structure 
over  the  Navier-Stokes  equations. 

The  Burnett  equations  are  derived  by  retaining  second 
order  terms  in  the  Chapman-Enskog  expansion  of  the  Boltz¬ 
mann  equation.  The  Navier-Stokes  equations  can  be  de¬ 
rived  in  a  similiar  fashion  by  retaining  only  first  order  terms 
in  the  Chapman-Enskog  expansion.  The  Burnett  equations 
differ  from  the  Navier-Stokes  equations  only  in  the  expres¬ 
sions  for  the  stress  tensor  and  the  heat  flux  vector.  The 
Burnett  stress  tensor  and  heat  flux  vector  are  the  Navier- 
Stokes  stress  tensor  and  heat  flux  vector  plus  several  addi¬ 
tional  terms.  The  Burnett  stress  tensor  and  heat  flux  vector 
in  one  spatial  dimension  for  a  hard  sphere  gas  are: 
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Due  to  the  complexity  of  the  Burnett  equations  a  sim¬ 
plified  constitutive  model  has  been  developed  based  on  an 
earlier  model  developed  by  Chapman,  et  al .'  This  model, 
hereafter  refered  to  as  the  nonlinear  constitutive  model,  at¬ 
tempts  to  capture  the  significant  features  of  the  Burnett 
constitutive  relations  with  only  two  additional  terms  to  the 
Navier-Stokes  constitutive  relations  as  shown  below. 
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)  The  additional  terms  in  the  stress  and  heat  flux  have  been 

chosen  such  that  the  model,  when  applied  to  shock  wave 
structure  in  argon,  gives  experimentally  measured  shock 
wave  reciprocal  thickness  and  approximately  the  same  sep¬ 
aration  between  the  temperature  and  density  profiles  calcu¬ 
lated  by  Fiscko  and  Chapman  for  argon  using  the  DSMC 
technique.  Fig.  7  shows  the  reciprocal  thicknesses  calcu¬ 
lated  using  the  nonlinear  constitutive  model  and  those  ex¬ 
perimentally  measured  by  Alsmeyer14 . 

The  nonlinear  constitutive  model  presented  above  dif¬ 
fers  from  the  the  model  presented  by  Chapman,  eial .}  only 
in  the  coefficient  of  the  velocity  gradient  squared  term  in  the 
stress,  which  has  been  changed  from  7  to  9 .5 .  This  change 
is  because  the  earlier  model  was  developed  to  give  the  cor¬ 
rect  shock  thicknesses  in  nitrogen  when  using  the  Landau- 
Teller  model  and  the  Parker  expression  for  Zr  adjusted  to 
fit  the  Lordi-Mates  computations.  Since  the  purpose  of  this 
paper  is  to  evaluate  rotational  models,  it  was  necessary  to 
develop  constitutive  relations  that  are  not  dependent  on  the 
rotational  model  employed.  Thus,  comparison  to  argon  re¬ 
ciprocal  thickness  was  chosen  to  calibrate  the  nonlinear 
constitutive  model.  It  is  then  assumed  that  the  nonlinear 
constitutive  model  will  be  as  effective  in  modeling  the  true 
stress  and  heat  flux  in  a  diatomic  gas  as  in  a  monatomic  gas . 

Application  of  Rotational  Noneauilibrium 
to  Shock  Wave  Structure 


As  mentioned  above  the  study  of  shock  wave  structure 
in  diatomic  gases  requires  the  inclusion  of  rotational  ther¬ 
mal  nonequilibrium.  A  method  which  “fully  couples”  these 
nonequilibrium  effects  into  the  equations  describing  the 
fluid  mechanics  was  chosen  in  this  study.  This  is  done  by 
including  a  rotational  energy  equation  into  the  continuum 
equations  of  motion.  In  order  to  do  this  the  relaxation  model 
must  be  transformed  from  the  stationary  forms  presented 
above  into  a  form  which  accounts  for  bulk  fluid  motion.  In 
addition  conduction  of  rotational  energy  must  be  included 
in  the  form  of  the  Fourier  Law  of  Heat  Conduction  in  both 
the  rotational  and  total  energy  equations.  Eq.(l)andEq.  (8) 
thus  become  Eq.  (14). 
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Note  that  the  ordinary  derivative  has  been  replaced  by  a 
substantial  derivative  to  account  for  convection  of  rota¬ 
tional  energy  due  to  fluid  motion.  Also  note  that  Eq.  (14a) 
is  no  longer  homogeneous,  but  now  has  a  thermal  source 
term.  Thus,  the  governing  equations  for  one  dimensional 
shock  structure  in  “conservation  law”  form  for  a  diatomic 
gas  are: 


(¥)’•  New  Model  b 
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The  method  used  to  solve  the  above  equation  set  is  sim- 
iliar  to  the  one  used  for  monatomic  gases  in  Ref.  (10)  which 
treats  the  Euler  terms  explicitly  and  the  viscous  terms  im¬ 
plicitly.  The  main  difference  here  is  the  non-homogeneous 
source  term  in  the  rotational  energy  equation.  The  equa¬ 
tion  set  for  the  monatomic  case  is  homogenous.  This  re¬ 
quires  a  slight  modification  to  the  numerical  technique  out¬ 
lined  in  Ref.  (10).  The  source  term  here  is  treated  in  an  im¬ 
plicit  manner  similiar  to  the  treatment  of  the  viscous  terms. 
Details  of  handling  source  terms  implicitly  are  outlined  in 
Refs.  (15,16). 

Shock  Structure  Results 


The  continuum  technique  described  above  was  used  to 
calculate  shock  wave  structure  for  four  combinations  of 
constitutive  model  and  rotational  model  at  eight  Mach  num¬ 
bers.  The  four  cases  considered  were  Burnett  constitu¬ 
tive  relations  (Eq.  12)  with  both  the  new  rotational  model 
and  the  Landau-Teller  model  (using  the  Lordi-Mates  values 
for  Zr)  and  the  nonlinear  model  constitutive  relations  (Eq. 
1 3)  with  the  same  two  rotational  models  considered  above. 
Shock  structures  were  calculated  at  Mach  1 .2, 1.5, 2, 3, 5, 6, 
8,  and  11.  In  all  cases  upstream  temperatures  were  chosen 
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10  be  in  equilibrium  and  at  300K,  and  upstream  pressure  was 
chosen  to  be  one  atm.  Shock  structure  density  and  temper¬ 
ature  profiles  were  normalized  as  shown  schematically  in 
Fig.  8.  The  normalization  values  pi—  p\  and  Ti  —T\  were 
calculated  analytically  with  the  assumption  that  the  down¬ 
stream  condition  was  vibrationally  frozen.  Table  4  summa¬ 
rizes  the  normalization  values  for  density  and  temperature 
jumps  used  in  this  work.  The  spatial  coordinate  was  nor¬ 
malized  by  the  upstream  mean  free  path  as  defined  below. 

Xs^Ijw,  (,6) 

Nitrogen  was  modeled  as  having  a  spherically  repulsive  po¬ 
tential.  This  gave  the  following  relation  between  viscosity 
and  temperature. 

H  =  iire/(=P-r  07) 

ire/ 

The  temperature  viscosity  exponent  was  chosen  to  be  0.72. 

Fig.  9  compares  the  rotational  and  translational  temper¬ 
ature  profiles  of  the  new  rotational  model  and  the  Landau- 
Teller  model  when  Burnett  constitutive  relations  are  used 
to  calculate  a  Mach  6  shock  structure.  Fig.  10  makes 
the  same  comparison  for  the  nonlinear  constitutive  rela¬ 
tions.  In  both  comparisons  the  new  rotational  model  gives  a 
lowerpeak  translational  temperature  than  the  Landau-Teller 
model.  However  the  profiles  in  both  comparisons  are  qual¬ 
itatively  very  similar.  This  suggests  that  the  Landau-Teller 
model  with  a  slightly  adjusted  value  of  ZR  may  reproduce 
the  results  acheived  with  the  new  rotational  model.  Figs.  1 1 
and  12  compare  the  computed  values  of  chock  reciprocal 
thickness  (as  defined  in  Fig.  8)  with  experimentally  mea¬ 
sured  values  from  Ref.  14.  Fiscko  and  Chapman10  show 
that  the  Burnett  constitutive  relations  yield  shock  recipro¬ 
cal  thicknesses  which  are  slightly  higher  than  the  experi¬ 
mental  data  for  argon;  therefore,  any  differences  between 
the  computed  and  experimental  shock  thicknesses  in  Fig.  1 1 
do  not  allow  a  rigorous  evaluation  of  the  rotational  model 
employed.  The  nonlinear  model  was  adjusted  to  give  the 
correct  shock  thickness  in  argon  over  the  range  of  Mach 
numbers  from  1  to  11  (Fig.  7);  therefore,  the  comparison  of 
computed  and  measured  reciprocal  thicknesses  in  Fig.  1 1 
should  provide  insight  as  to  the  validity  of  the  rotational 
models  employed. 

The  new  rotational  model  accurately  predicts  shock  re¬ 
ciprocal  thickness  at  Mach  numbers  below  six;  however, 
above  Mach  6  the  new  rotational  model  over  predicts  shock 
reciprocal  thickness.  The  Landau-Teller  model  using  the 
Lordi-Mates  values  for  Z&  agrees  with  experiment  to  only 
Mach  3  and  under  predicts  reciprocal  thickness  thereafter. 


The  deficiency  in  the  n.w  model  *s  thought  to  be  in  the 
measured  transition  rates  that  were  used  to  develop  the 
model.  Rahn, era/..1112,13  measured  the  transition  rates  up 
to  2000K;  therefore,  the  curve  fits  suggested  for  the  tran¬ 
sition  rates  may  be  in  error  for  temperatures  higher  than 
2000K.  This  would  explain  the  inaccuraccy  in  the  new  ro¬ 
tational  model  at  Mach  numbers  higher  than  six  where  the 
downstream  temperature  is  in  excess  of  2000K. 

Finally  in  view  of  the  inaccuracies  of  both  rotational 
models  at  high  Mach  numbers,  the  expression  for  ZR  in 
the  Landau-Teller  equation  was  modified  in  an  attempt  to 
give  correct  shock  reciprocal  thicknesses.  The  value  of  ZR 
was  adjusted  to  18.0  while  T *  was  maintained  at  91. 5K. 
The  results  of  this  simple  modification  are  also  presented 
in  Fig.  12  and  agree  well  with  the  experimental  data  at  all 
Mach  numbers  up  to  1 1 . 

Conclusion 

In  conclusion,  the  Landau-Teller  equation  is  presently 
being  used  to  model  rotational  relaxation  in  hypersonic 
CFD,  even  though  the  assumptions  used  to  derive  the 
Landau-Teller  equation  are  probably  violated  in  hypersonic 
flow  fields.  However,  it  has  been  found  that  by  adjusting  pa¬ 
rameters  in  Parker’s  analytical  expression  for  the  rotational 
collision  number  ZR,  the  Landau-Teller  model  along  with  a 
nonlinear  constitutive  model  will  accurately  reproduce  ex¬ 
perimental  shock  reciprocal  thicknesses  for  nitrogen.  Thus, 
while  neither  Landau-Teller  relaxation,  nor  ZR  =  18  .0 
in  Parker’s  expression  may  be  individually  accurate  for  ni¬ 
trogen,  the  combination  of  the  two  yields  accurate  shock 
thicknesses  up  to  Mach  11,  the  limit  of  the  experimental 
data.  Due  to  a  lack  of  experimental  shock  wave  tempera¬ 
ture  profiles  in  nitrogen,  it  is  unclear  how  well  this  method 
reproduces  the  detailed  structure  of  the  nitrogen  shock. 

In  an  effort  to  more  realistically  capture  the  true  physics 
of  rotational  thermal  nonequilibrium,  an  alternative  to  the 
Landau-Teller  equation  which  is  based  on  the  more  rigorous 
master  equation  has  been  developed.  This  alternative,  in 
contrast  to  Landau-Teller,  yields  a  rotational  collision  num¬ 
ber  (Zg)  dependent  on  Tr /Tt  as  well  as  Tt ,  and  a  time  rate  of 
rotational  energy  change  that  is  proportional  to  some  power 
(other  than  one)  of  the  energy  difference,  e\q(Tt)  -eT.  This 
new  model  has  accurately  reproduced  experimental  shock 
reciprocal  thicknesses  in  nitrogen  up  to  Mach  6  and  slightly 
over  predicts  shock  reciprocal  thickness  above  Mach  6.  It  is 
believed  that  the  deficiency  in  the  new  mode!  is  the  uncer- 
tainity  at  the  high  temperatures  found  in  these  high  Mach 
number  shocks  of  the  transition  rates-extrapolatcd  above 
2000K-upon  which  the  model  is  based. 
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Finally,  n  should  be  noted  that  while  the  rotational  tran¬ 
sition  rates  of  Rahn,  et  a/.,111213  have  an  applicable  tem¬ 
perature  ranee  of  under  2000K,  the  applicability  of  this 
approach  ot  so  limited.  As  these  rates  become  better 
known  though  the  efforts  of  experimentalists  and  compu¬ 
tational  chemists,  the  improved  values  can  be  similarly  used 
to  arrive  at  even  more  refined  models  for  rotational  thermal 
nonequilibrium. 
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Table  1.  Values  used  in  Parker’s  Model. 


T'  (K) 

2? 

Parker 

80 

15.7 

Lordi/Mates 

91.5 

23.0 

Table  2.  Values  used  in  empirical  fit  of  /f, 


a  (Pa  sec)-1 

0 

6 

m  |  T0  (K) 

6805.1 

1.67 

1.21 

0.1487  j  295 

Table  3.  Values  used  in  fits  for  Zr  and  n. 


i 

1 

2 

3 

4 

5 

I  6 

At 

6.16574 

3.9166 

7.919 

2.1858 

5.0616 

I  0.6035 

Bi 

Hffig&QEBU 

1.8 

1250.0 

500.0 

I 

C, 

1  1.1282934 

-0.04920 

0.762538 

-0.15419 

Table  4.  Values  used  to  normalize  shock  wave  profiles. 


Mach 

Pi  (kg/m-1 j 

Pi  -  Pi  (kg/m-1) 

Tz  (K) 

Tz  -  fi  (K) 

1.2 

1.526 

0.389 

338.4 

38.44 

1.5 

2.119 

0.981 

396.1 

96.07 

2.0 

3.035 

1.897 

506.3 

206.3 

3.0 

4.390 

3.252 

803.7 

503.7 

5.0 

5.691 

4.553 

1740.0 

1440.0 

6.0 

5.997 

4.858 

2382.2 

2082.2 

8.0 

6.334 

5.196 

4016.0 

3716.0 

11.0 

6.558 

5.42C 

7341.3 

7041.3 

Model  Exponent  (n)  Temperature  (Kelvin) 
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Fig.  1  Rotational  collision  number  vs.  temperature. 


Fig.  4  New  rotational  model  collision  number  Zr,  Tt  <  Tt. 
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Fig.  2  Adiabatic  box  temperature  relaxation  time  histories.  Fig.  5  New  rotational  model  exponent  n,  Tr  >  Tt. 


Fig.  6  New  rotational  model  collision  number  Zr.T,  >  Tt. 


Fig.  3  New  rotational  model  exponent  n,  7V  <Tt. 
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Fig.  7  Shock  reciprocal  thickness  in  argon.  Comparison  of 
experiment  and  results  using  nonlinear  constiutive  model. 


Fig.  8  Schematic  of  shock  structure  parameters  and  nor¬ 
malizations. 


Fig.  9  Mach  6  temperature  profiles  using  Burnett  constitu¬ 
tive  relations.  Nitrogen,  w  =  0 .72 . 


Fig.  10  Mach  6  temperature  profiles  using  nonlinear  model 
constitutive  relations.  Nitrogen,  ui  =  0 .72 . 


Fig.  11  Shock  reciprocal  thickness  -  comparison  of  Burnett 
results  and  experiment.  Nitrogen,  w  =  0 .72 . 


Fig.  12  Shock  reciprocal  thickness  -  comparison  of  nonlin¬ 
ear  model  results  and  experiment.  Nitrogen. «  =  0 .72 . 


